package net.yegong.matrix;

import org.netlib.lapack.Sgetrf;
import org.netlib.util.intW;

/**
 * @author cooper
 */
public class LUDecomposition {

	public LUDecomposition(Matrix mat) {
		this.mat = mat;
	}

	public void calc() {
		int m = mat.nRow;
		int n = mat.nCol;
		int min = Math.min(m, n);
		int[] ipiv = new int[min];
		float[] a = mat.getColumnMajorArrayCopy();
		intW info = new intW(0);
		Sgetrf.sgetrf(m, n, a, 0, m, ipiv, 0, info);
		if (m == n) {
			l = new SquareMatrix(a, m, Matrix.FORTRAN);
		}
		else {
			l = new Matrix(a, m, n, Matrix.FORTRAN);
		}
		u = l.clone();
		l.reserveLowerTriangle();
		u.reserveUpperTriangle();
		for (int i = 0; i < min; ++i) {
			l.set(i, i, 1);
		}
		if (m > n) {// slice u to Square
			u = u.rowSlice(0, n);
		}
		if (m < n) {// slice l to Square
			l = l.columnSlice(0, m);
		}
		int[] map = new int[m];
		for (int i = 0; i < m; ++i) {
			map[i] = i;
		}
		for (int i = 0; i < min; ++i) {
			--ipiv[i];
			int tmp = map[i];
			map[i] = map[ipiv[i]];
			map[ipiv[i]] = tmp;
		}
		p = new PermutationMatrix(map, PermutationMatrix.FROM_MAP, PermutationMatrix.ROWS);
	}

	public float det() {
		if (u == null)
			calc();
		return u.getDiag().det() * p.det();
	}

	public Matrix getL() {
		if (l == null)
			calc();
		return l;
	}

	public Matrix getU() {
		if (u == null)
			calc();
		return u;
	}

	public PermutationMatrix getP() {
		if (p == null)
			calc();
		return p;
	}

	protected Matrix mat;
	protected Matrix l;
	protected Matrix u;
	protected PermutationMatrix p;
}
